{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/NeuromatchAcademy/course-content/blob/master/tutorials/W3D1_BayesianDecisions/W3D1_Tutorial3.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Neuromatch Academy: Week 3, Day 1, Bonus Tutorial\n",
    "# Fitting to data\n",
    "\n",
    "__Content creators:__ Vincent Valton, Konrad Kording\n",
    "\n",
    "__Content reviewers:__ Matt Krause, Jesse Livezey, Karolina Stosio, Saeed Salehi, Michael Waskom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Tutorial objectives\n",
    "  \n",
    "In the first two tutorials, we learned about Bayesian models and decisions more intuitively, using demos. In this notebook, we will dive into using math and code to fit Bayesian models to data. \n",
    "\n",
    "We'll have a look at computing all the necessary steps to perform model inversion (estimate the model parameters such as $p_{common}$ that generated data similar to that of a participant). We will describe all the steps of the generative model first, and in the last exercise we will use all these steps to estimate the parameter $p_{common}$ of a single participant using simulated data.   \n",
    "\n",
    "The generative model will be a Bayesian model we saw in Tutorial 2: a mixture of Gaussian prior  and a Gaussian likelihood.\n",
    "Steps:\n",
    "\n",
    "* First, we'll create the prior, likelihood, posterior, etc in a form that will make it easier for us to visualise what is being computed and estimated at each step of the generative model: \n",
    "  1. Creating a mixture of Gaussian prior for multiple possible stimulus inputs\n",
    "  2. Generating the likelihood for multiple possible stimulus inputs\n",
    "  3. Estimating our posterior as a function of the stimulus input\n",
    "  4. Estimating a participant response given the posterior\n",
    "  \n",
    "* Next, we'll perform the model inversion/fitting:\n",
    "  5. Create an distribution for the input as a function of possible inputs\n",
    "  6. Marginalization\n",
    "  7. Generate some data using the generative model provided\n",
    "  8. Perform model inversion (model fitting) using the generated data and see if you recover the orignal parameters.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Setup\n",
    "\n",
    "Please execute the cell below to initialize the notebook environment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:51.013866Z",
     "iopub.status.busy": "2021-06-02T13:21:51.013268Z",
     "iopub.status.idle": "2021-06-02T13:21:54.019087Z",
     "shell.execute_reply": "2021-06-02T13:21:54.018504Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "from scipy.optimize import minimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:54.025116Z",
     "iopub.status.busy": "2021-06-02T13:21:54.023450Z",
     "iopub.status.idle": "2021-06-02T13:21:54.217843Z",
     "shell.execute_reply": "2021-06-02T13:21:54.217321Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Figure Settings\n",
    "import ipywidgets as widgets\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "plt.style.use(\"https://raw.githubusercontent.com/NeuromatchAcademy/course-content/NMA2020/nma.mplstyle\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:54.231001Z",
     "iopub.status.busy": "2021-06-02T13:21:54.225546Z",
     "iopub.status.idle": "2021-06-02T13:21:54.233228Z",
     "shell.execute_reply": "2021-06-02T13:21:54.233589Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title Helper Functions\n",
    "\n",
    "def my_gaussian(x_points, mu, sigma):\n",
    "    \"\"\"\n",
    "    Returns a Gaussian estimated at points `x_points`, with parameters: `mu` and `sigma`\n",
    "\n",
    "    Args :\n",
    "      x_points (numpy arrays of floats)- points at which the gaussian is evaluated\n",
    "      mu (scalar) - mean of the Gaussian\n",
    "      sigma (scalar) - std of the gaussian\n",
    "\n",
    "    Returns:\n",
    "      Gaussian evaluated at `x`\n",
    "    \"\"\"\n",
    "    p = np.exp(-(x_points-mu)**2/(2*sigma**2))\n",
    "    return p / sum(p)\n",
    "\n",
    "def moments_myfunc(x_points, function):\n",
    "  \"\"\"\n",
    "  DO NOT EDIT THIS FUNCTION !!!\n",
    "\n",
    "  Returns the mean, median and mode of an arbitrary function\n",
    "\n",
    "  Args :\n",
    "    x_points (numpy array of floats) - x-axis values\n",
    "    function (numpy array of floats) - y-axis values of the function evaluated at `x_points`\n",
    "\n",
    "  Returns:\n",
    "    (tuple of 3 scalars): mean, median, mode\n",
    "  \"\"\"\n",
    "\n",
    "  # Calc mode of arbitrary function\n",
    "  mode = x_points[np.argmax(function)]\n",
    "\n",
    "  # Calc mean of arbitrary function\n",
    "  mean = np.sum(x_points * function)\n",
    "\n",
    "  # Calc median of arbitrary function\n",
    "  cdf_function = np.zeros_like(x_points)\n",
    "  accumulator = 0\n",
    "  for i in np.arange(x_points.shape[0]):\n",
    "    accumulator = accumulator + function[i]\n",
    "    cdf_function[i] = accumulator\n",
    "  idx = np.argmin(np.abs(cdf_function - 0.5))\n",
    "  median = x_points[idx]\n",
    "\n",
    "  return mean, median, mode\n",
    "\n",
    "def plot_myarray(array, xlabel, ylabel, title):\n",
    "  \"\"\" Plot an array with labels.\n",
    "\n",
    "  Args :\n",
    "    array (numpy array of floats)\n",
    "    xlabel (string) - label of x-axis\n",
    "    ylabel (string) - label of y-axis\n",
    "    title  (string) - title of plot\n",
    "\n",
    "  Returns:\n",
    "    None\n",
    "  \"\"\"\n",
    "  fig = plt.figure()\n",
    "  ax = fig.add_subplot(111)\n",
    "  colormap = ax.imshow(array, extent=[-10, 10, 8, -8])\n",
    "  cbar = plt.colorbar(colormap, ax=ax)\n",
    "  cbar.set_label('probability')\n",
    "  ax.invert_yaxis()\n",
    "  ax.set_xlabel(xlabel)\n",
    "  ax.set_title(title)\n",
    "  ax.set_ylabel(ylabel)\n",
    "  ax.set_aspect('auto')\n",
    "  return None\n",
    "\n",
    "def plot_my_bayes_model(model) -> None:\n",
    "  \"\"\"Pretty-print a simple Bayes Model (ex 7), defined as a function:\n",
    "\n",
    "  Args:\n",
    "    - model: function that takes a single parameter value and returns\n",
    "             the negative log-likelihood of the model, given that parameter\n",
    "  Returns:\n",
    "    None, draws plot\n",
    "    \"\"\"\n",
    "  x = np.arange(-10,10,0.07)\n",
    "\n",
    "  # Plot neg-LogLikelihood for different values of alpha\n",
    "  alpha_tries = np.arange(0.01, 0.3, 0.01)\n",
    "  nll = np.zeros_like(alpha_tries)\n",
    "  for i_try in np.arange(alpha_tries.shape[0]):\n",
    "    nll[i_try] = model(np.array([alpha_tries[i_try]]))\n",
    "\n",
    "  plt.figure()\n",
    "  plt.plot(alpha_tries, nll)\n",
    "  plt.xlabel('p_independent value')\n",
    "  plt.ylabel('negative log-likelihood')\n",
    "\n",
    "  # Mark minima\n",
    "  ix = np.argmin(nll)\n",
    "  plt.scatter(alpha_tries[ix], nll[ix], c='r', s=144)\n",
    "\n",
    "  #plt.axvline(alpha_tries[np.argmin(nll)])\n",
    "  plt.title('Sample Output')\n",
    "  plt.show()\n",
    "\n",
    "  return None\n",
    "\n",
    "\n",
    "def plot_simulated_behavior(true_stim, behaviour):\n",
    "  fig = plt.figure(figsize=(7, 7))\n",
    "  ax = fig.add_subplot(1,1,1)\n",
    "  ax.set_facecolor('xkcd:light grey')\n",
    "  plt.plot(true_stim, true_stim - behaviour, '-k', linewidth=2, label='data')\n",
    "  plt.axvline(0, ls='dashed', color='grey')\n",
    "  plt.axhline(0, ls='dashed', color='grey')\n",
    "  plt.legend()\n",
    "  plt.xlabel('Position of true visual stimulus (cm)')\n",
    "  plt.ylabel('Participant deviation from true stimulus (cm)')\n",
    "  plt.title('Participant behavior')\n",
    "  plt.show()\n",
    "\n",
    "  return None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Introduction\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:54.240506Z",
     "iopub.status.busy": "2021-06-02T13:21:54.240023Z",
     "iopub.status.idle": "2021-06-02T13:21:54.333035Z",
     "shell.execute_reply": "2021-06-02T13:21:54.333425Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Video 1: Intro\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id='YSKDhnbjKmA', width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtube.com/watch?v=\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "![Generative model](https://github.com/vincentvalton/figures_NMA_W2D1_T3/blob/master/Drawing%20Generative%20Model%20W2T3.png?raw=true)\n",
    "\n",
    "Here is a graphical representation of the generative model:\n",
    "\n",
    "  1. We present a stimulus $x$ to participants. \n",
    "  2. The brain encodes this true stimulus $x$ noisily (this is the brain's representation of the true visual stimulus: $p(\\tilde x|x)$.\n",
    "  3. The brain then combine this brain encoded stimulus (likelihood: $p(\\tilde x|x)$) with prior information (the prior: $p(x)$) to make up the brain's estimated position of the true visual stimulus, the posterior: $p(x|\\tilde x)$. \n",
    "  3. This brain's estimated stimulus position: $p(x|\\tilde x)$, is then used to make a response:  $\\hat x$, which is the participant's noisy estimate of the stimulus position (the participant's percept). \n",
    "  \n",
    "Typically the response $\\hat x$ also includes some motor noise (noise due to the hand/arm move being not 100% accurate), but we'll ignore it in this tutorial and assume there is no motor noise.\n",
    "\n",
    "\n",
    "\n",
    "We will use the same experimental setup as in [tutorial 2](https://colab.research.google.com/drive/15pbgrfGjSKbUQoX51RdcNe3UXb4R5RRx#scrollTo=tF5caxVGYURh) but with slightly different probabilities. This time, participants are told that they need to estimate the sound location of a puppet that is hidden behind a curtain. The participants are told to use auditory information and are also informed that the sound could come from 2 possible causes: a common cause (95% of the time it comes from the puppet hidden behind the curtain at position 0), or an independent cause (5% of the time the sound comes from loud-speakers at more distant locations)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Section 1: Likelihood array\n",
    "    \n",
    "First, we want to create a likelihood, but for the sake of visualization (and to consider all possible brain encodings) we will create multiple likelihoods $f(x)=p(\\tilde x|x)$ (one for each potential encoded stimulus: $\\tilde x$). We will then be able to visualize the likelihood as a function of hypothesized true stimulus positions: $x$ on the x-axis and encoded position $\\tilde x$ on the y-axis.\n",
    "\n",
    "\n",
    "  Using the equation for the `my_gaussian` and the values in `hypothetical_stim`:\n",
    "* Create a Gaussian likelihood with mean varying from `hypothetical_stim`, keeping $\\sigma_{likelihood}$ constant at 1.\n",
    "* Each likelihood will have a different mean and thus a different row-likelihood of your 2D array, such that you end up with a likelihood array made up of 1,000 row-Gaussians with different means. (_Hint_: `np.tile` won't work here. You may need a for-loop).\n",
    "* Plot the array using the function `plot_myarray()` already pre-written and commented-out in your script"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Exercise 1. Implement the auditory likelihood as a function of true stimulus position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "code",
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:54.339528Z",
     "iopub.status.busy": "2021-06-02T13:21:54.338470Z",
     "iopub.status.idle": "2021-06-02T13:21:54.340079Z",
     "shell.execute_reply": "2021-06-02T13:21:54.340463Z"
    }
   },
   "outputs": [],
   "source": [
    "x = np.arange(-10, 10, 0.1)\n",
    "hypothetical_stim = np.linspace(-8, 8, 1000)\n",
    "\n",
    "def compute_likelihood_array(x_points, stim_array, sigma=1.):\n",
    "\n",
    "    # initializing likelihood_array\n",
    "    likelihood_array = np.zeros((len(stim_array), len(x_points)))\n",
    "\n",
    "    # looping over stimulus array\n",
    "    for i in range(len(stim_array)):\n",
    "        ########################################################################\n",
    "        ## Insert your code here to:\n",
    "        ##      - Generate a likelihood array using `my_gaussian` function,\n",
    "        ##        with std=1, and varying the mean using `stim_array` values.\n",
    "        ## remove the raise below to test your function\n",
    "        raise NotImplementedError(\"You need to complete the function!\")\n",
    "        ########################################################################\n",
    "        likelihood_array[i, :] = ...\n",
    "\n",
    "    return likelihood_array\n",
    "\n",
    "# Uncomment following lines to test your code\n",
    "# likelihood_array = compute_likelihood_array(x, hypothetical_stim)\n",
    "# plot_myarray(likelihood_array,\n",
    "#             '$x$ : Potential true stimulus $x$',\n",
    "#             'Possible brain encoding $\\~x$',\n",
    "#             'Likelihood as a function of $\\~x$ : $p(\\~x | x)$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:54.348748Z",
     "iopub.status.busy": "2021-06-02T13:21:54.345280Z",
     "iopub.status.idle": "2021-06-02T13:21:55.254562Z",
     "shell.execute_reply": "2021-06-02T13:21:55.254978Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "x = np.arange(-10, 10, 0.1)\n",
    "hypothetical_stim = np.linspace(-8, 8, 1000)\n",
    "\n",
    "def compute_likelihood_array(x_points, stim_array, sigma=1.):\n",
    "\n",
    "    # initializing likelihood_array\n",
    "    likelihood_array = np.zeros((len(stim_array), len(x_points)))\n",
    "\n",
    "    # looping over stimulus array\n",
    "    for i in range(len(stim_array)):\n",
    "        likelihood_array[i, :] = my_gaussian(x_points, stim_array[i], sigma)\n",
    "\n",
    "    return likelihood_array\n",
    "\n",
    "likelihood_array = compute_likelihood_array(x, hypothetical_stim)\n",
    "\n",
    "with plt.xkcd():\n",
    "  plot_myarray(likelihood_array,\n",
    "              '$x$ : Potential true stimulus $x$',\n",
    "              'Possible brain encoding $\\~x$',\n",
    "              'Likelihood as a function of $\\~x$ : $p(\\~x | x)$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Section 2: Causal mixture of Gaussian prior\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:55.259896Z",
     "iopub.status.busy": "2021-06-02T13:21:55.259388Z",
     "iopub.status.idle": "2021-06-02T13:21:55.346948Z",
     "shell.execute_reply": "2021-06-02T13:21:55.346517Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Video 2: Prior array\n",
    "video = YouTubeVideo(id='F0IYpUicXu4', width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtube.com/watch?v=\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "As in Tutorial 2, we want to create a prior that will describe the participants' prior knowledge that, 95% of the time sounds come from a common position around the puppet, while during the remaining 5% of the time, they arise from another independent position. We will embody this information into a prior using a mixture of Gaussians. For visualization reasons, we will create a prior that has the same shape (form) as the likelihood array we created in the previous exercise. That is, we want to create a mixture of Gaussian prior as a function the the brain encoded stimulus $\\tilde x$. Since the prior does not change as a function of $\\tilde x$ it will be identical for each row of the prior 2D array. \n",
    "\n",
    "Using the equation for the Gaussian `my_gaussian`:\n",
    "* Generate a Gaussian $Common$ with mean 0 and standard deviation 0.5\n",
    "* Generate another Gaussian $Independent$ with mean 0 and standard deviation 10\n",
    "* Combine the two Gaussians (Common + Independent) to make a new prior by mixing the two Gaussians with mixing parameter $p_{independent}$ = 0.05. Make it such that the peakier Gaussian has 95% of the weight (don't forget to normalize afterwards)\n",
    "* This will be the first row of your prior 2D array\n",
    "* Now repeat this for varying brain encodings $\\tilde x$. Since the prior does not depend on $\\tilde x$ you can just repeat the prior for each $\\tilde x$ (hint: use np.tile) that row prior to make an array of 1,000 (i.e. `hypothetical_stim.shape[0]`)  row-priors.\n",
    "* Plot the matrix using the function `plot_myarray()` already pre-written and commented-out in your script"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2: Implement the prior array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "code",
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:55.353160Z",
     "iopub.status.busy": "2021-06-02T13:21:55.352029Z",
     "iopub.status.idle": "2021-06-02T13:21:55.353716Z",
     "shell.execute_reply": "2021-06-02T13:21:55.354103Z"
    }
   },
   "outputs": [],
   "source": [
    "x = np.arange(-10, 10, 0.1)\n",
    "\n",
    "def calculate_prior_array(x_points, stim_array, p_indep,\n",
    "                          prior_mean_common=.0, prior_sigma_common=.5,\n",
    "                          prior_mean_indep=.0, prior_sigma_indep=10):\n",
    "    \"\"\"\n",
    "        'common' stands for common\n",
    "        'indep' stands for independent\n",
    "    \"\"\"\n",
    "\n",
    "    prior_common = my_gaussian(x_points, prior_mean_common, prior_sigma_common)\n",
    "    prior_indep = my_gaussian(x_points, prior_mean_indep, prior_sigma_indep)\n",
    "\n",
    "    ############################################################################\n",
    "    ## Insert your code here to:\n",
    "    ##      - Create a mixture of gaussian priors from 'prior_common'\n",
    "    ##        and 'prior_indep' with mixing parameter 'p_indep'\n",
    "    ##      - normalize\n",
    "    ##      - repeat the prior array and reshape it to make a 2D array\n",
    "    ##        of 1000 rows of priors (Hint: use np.tile() and np.reshape())\n",
    "    ## remove the raise below to test your function\n",
    "    raise NotImplementedError(\"You need to complete the function!\")\n",
    "    ############################################################################\n",
    "\n",
    "    prior_mixed = ...\n",
    "    prior_mixed /= ...  # normalize\n",
    "\n",
    "    prior_array = np.tile(...).reshape(...)\n",
    "    return prior_array\n",
    "\n",
    "p_independent=.05\n",
    "# Uncomment following lines, once the task is complete.\n",
    "# prior_array = calculate_prior_array(x, hypothetical_stim, p_independent)\n",
    "# plot_myarray(prior_array,\n",
    "#              'Hypothesized position $x$', 'Brain encoded position $\\~x$',\n",
    "#              'Prior as a fcn of $\\~x$ : $p(x|\\~x)$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:55.377197Z",
     "iopub.status.busy": "2021-06-02T13:21:55.376705Z",
     "iopub.status.idle": "2021-06-02T13:21:55.974053Z",
     "shell.execute_reply": "2021-06-02T13:21:55.974430Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "x = np.arange(-10, 10, 0.1)\n",
    "\n",
    "def calculate_prior_array(x_points, stim_array, p_indep,\n",
    "                          prior_mean_common=.0, prior_sigma_common=.5,\n",
    "                          prior_mean_indep=.0, prior_sigma_indep=10):\n",
    "    \"\"\"\n",
    "        'common' stands for common\n",
    "        'indep' stands for independent\n",
    "    \"\"\"\n",
    "\n",
    "    prior_common = my_gaussian(x_points, prior_mean_common, prior_sigma_common)\n",
    "    prior_indep = my_gaussian(x_points, prior_mean_indep, prior_sigma_indep)\n",
    "\n",
    "    prior_mixed = (1 - p_indep) * prior_common + (p_indep * prior_indep)\n",
    "    prior_mixed /= np.sum(prior_mixed)  # normalize\n",
    "\n",
    "    prior_array = np.tile(prior_mixed, len(stim_array)).reshape(len(stim_array), -1)\n",
    "    return prior_array\n",
    "\n",
    "p_independent=.05\n",
    "prior_array = calculate_prior_array(x, hypothetical_stim, p_independent)\n",
    "\n",
    "with plt.xkcd():\n",
    "  plot_myarray(prior_array,\n",
    "               'Hypothesized position $x$', 'Brain encoded position $\\~x$',\n",
    "               'Prior as a fcn of $\\~x$ : $p(x|\\~x)$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Section 3: Bayes rule and Posterior array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:55.979508Z",
     "iopub.status.busy": "2021-06-02T13:21:55.978990Z",
     "iopub.status.idle": "2021-06-02T13:21:56.068634Z",
     "shell.execute_reply": "2021-06-02T13:21:56.069030Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Video 3: Posterior array\n",
    "video = YouTubeVideo(id='HpOzXZUKFJc', width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtube.com/watch?v=\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now want to calcualte the posterior using *Bayes Rule*. Since we have already created a likelihood and a prior for each brain encoded position $\\tilde x$, all we need to do is to multiply them row-wise. That is, each row of the posterior array will be the posterior resulting from the multiplication of the prior and likelihood of the same equivalent row.\n",
    "\n",
    "Mathematically:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "    Posterior\\left[i, :\\right] \\propto Likelihood\\left[i, :\\right] \\odot Prior\\left[i, :\\right]\n",
    "\\end{eqnarray}\n",
    "\n",
    "where $\\odot$ represents the [Hadamard Product](https://en.wikipedia.org/wiki/Hadamard_product_(matrices)) (i.e., elementwise multiplication) of the corresponding prior and likelihood row vectors `i` from each matrix.\n",
    "\n",
    "Follow these steps to build the posterior as a function of the brain encoded stimulus $\\tilde x$:\n",
    "* For each row of the prior and likelihood (i.e. each possible brain encoding $\\tilde x$), fill in the posterior matrix so that every row of the posterior array represents the posterior density for a different brain encode  $\\tilde x$.\n",
    "* Plot the array using the function `plot_myarray()` already pre-written and commented-out in your script\n",
    "\n",
    "Optional:\n",
    "* Do you need to operate on one element--or even one row--at a time? NumPy operations can often process an entire matrix in a single \"vectorized\" operation. This approach is often much faster and much easier to read than an element-by-element calculation.  Try to write a vectorized version that calculates the posterior without using any for-loops. _Hint_: look at `np.sum` and its keyword arguments."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3: Calculate the posterior as a function of the hypothetical stimulus x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "code",
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:56.073238Z",
     "iopub.status.busy": "2021-06-02T13:21:56.072754Z",
     "iopub.status.idle": "2021-06-02T13:21:56.075509Z",
     "shell.execute_reply": "2021-06-02T13:21:56.075869Z"
    }
   },
   "outputs": [],
   "source": [
    "def calculate_posterior_array(prior_array, likelihood_array):\n",
    "\n",
    "    ############################################################################\n",
    "    ## Insert your code here to:\n",
    "    ##      - calculate the 'posterior_array' from the given\n",
    "    ##        'prior_array', 'likelihood_array'\n",
    "    ##      - normalize\n",
    "    ## remove the raise below to test your function\n",
    "    raise NotImplementedError(\"You need to complete the function!\")\n",
    "    ############################################################################\n",
    "    posterior_array = ...\n",
    "    posterior_array /= ...  # normalize each row separately\n",
    "\n",
    "    return posterior_array\n",
    "\n",
    "# Uncomment following lines, once the task is complete.\n",
    "# posterior_array = calculate_posterior_array(prior_array, likelihood_array)\n",
    "# plot_myarray(posterior_array,\n",
    "#              'Hypothesized Position $x$',\n",
    "#              'Brain encoded Stimulus $\\~x$',\n",
    "#              'Posterior as a fcn of $\\~x$ : $p(x | \\~x)$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:56.079860Z",
     "iopub.status.busy": "2021-06-02T13:21:56.079377Z",
     "iopub.status.idle": "2021-06-02T13:21:56.772373Z",
     "shell.execute_reply": "2021-06-02T13:21:56.771954Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "def calculate_posterior_array(prior_array, likelihood_array):\n",
    "\n",
    "    posterior_array = prior_array * likelihood_array\n",
    "    posterior_array /= posterior_array.sum(axis=1, keepdims=True)  # normalize each row separately\n",
    "\n",
    "    return posterior_array\n",
    "\n",
    "posterior_array = calculate_posterior_array(prior_array, likelihood_array)\n",
    "with plt.xkcd():\n",
    "  plot_myarray(posterior_array,\n",
    "               'Hypothesized Position $x$',\n",
    "               'Brain encoded Stimulus $\\~x$',\n",
    "               'Posterior as a fcn of $\\~x$ : $p(x | \\~x)$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Section 4: Estimating the position $\\hat x$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:56.777737Z",
     "iopub.status.busy": "2021-06-02T13:21:56.777261Z",
     "iopub.status.idle": "2021-06-02T13:21:56.895385Z",
     "shell.execute_reply": "2021-06-02T13:21:56.895838Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Video 4: Binary decision matrix\n",
    "video = YouTubeVideo(id='gy3GmlssHgQ', width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtube.com/watch?v=\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have a posterior distribution (for each possible brain encoding $\\tilde x$)that represents the brain's estimated stimulus position: $p(x|\\tilde x)$, we want to make an estimate (response) of the sound location $\\hat x$ using the posterior distribution. This would represent the subject's estimate if their (for us as experimentalist unobservable) brain encoding took on each possible value. \n",
    "\n",
    "This effectively encodes the *decision* that a participant would make for a given brain encoding $\\tilde x$. In this exercise, we make the assumptions that participants take the mean of the posterior (decision rule) as a response estimate for the sound location (use the function `moments_myfunc()` provided to calculate the mean of the posterior).\n",
    "\n",
    "Using this knowledge, we will now represent $\\hat x$ as a function of the encoded stimulus $\\tilde x$. This will result in a 2D binary decision array. To do so, we will scan the posterior matrix (i.e. row-wise), and set the array cell value to 1 at the mean of the row-wise posterior.\n",
    "\n",
    "**Suggestions**\n",
    "* For each brain encoding $\\tilde x$ (row of the posterior array), calculate the mean of the posterior, and set the corresponding cell of the binary decision array to 1. (e.g., if the mean of the posterior is at position 0, then set the cell with x_column == 0 to 1).\n",
    "* Plot the matrix using the function `plot_myarray()` already pre-written and commented-out in your script"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 4: Calculate the estimated response as a function of the hypothetical stimulus x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "code",
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:56.901541Z",
     "iopub.status.busy": "2021-06-02T13:21:56.900553Z",
     "iopub.status.idle": "2021-06-02T13:21:56.902159Z",
     "shell.execute_reply": "2021-06-02T13:21:56.902552Z"
    }
   },
   "outputs": [],
   "source": [
    "def calculate_binary_decision_array(x_points, posterior_array):\n",
    "\n",
    "    binary_decision_array = np.zeros_like(posterior_array)\n",
    "\n",
    "    for i in range(len(posterior_array)):\n",
    "\n",
    "        ########################################################################\n",
    "        ## Insert your code here to:\n",
    "        ##      - For each hypothetical stimulus x (row of posterior),\n",
    "        ##        calculate the mean of the posterior using the povided function\n",
    "        ##        `moments_myfunc()`, and set the corresponding cell of the\n",
    "        ##        Binary Decision array to 1.\n",
    "        ##        Hint: you can run 'help(moments_myfunc)' to see the docstring\n",
    "        ## remove the raise below to test your function\n",
    "        raise NotImplementedError(\"You need to complete the function!\")\n",
    "        ########################################################################\n",
    "        # calculate mean of posterior using 'moments_myfunc'\n",
    "        mean, _, _ = ...\n",
    "        # find the postion of mean in x_points (closest position)\n",
    "        idx = ...\n",
    "        binary_decision_array[i, idx] = 1\n",
    "\n",
    "    return binary_decision_array\n",
    "\n",
    "# Uncomment following lines, once the task is complete.\n",
    "# binary_decision_array = calculate_binary_decision_array(x, posterior_array)\n",
    "# plot_myarray(binary_decision_array,\n",
    "#              'Chosen position $\\hat x$', 'Brain-encoded Stimulus $\\~ x$',\n",
    "#              'Sample Binary Decision Array')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:56.907822Z",
     "iopub.status.busy": "2021-06-02T13:21:56.907061Z",
     "iopub.status.idle": "2021-06-02T13:21:57.609814Z",
     "shell.execute_reply": "2021-06-02T13:21:57.609395Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "def calculate_binary_decision_array(x_points, posterior_array):\n",
    "\n",
    "    binary_decision_array = np.zeros_like(posterior_array)\n",
    "\n",
    "    for i in range(len(posterior_array)):\n",
    "        # calculate mean of posterior using 'moments_myfunc'\n",
    "        mean, _, _ = moments_myfunc(x_points, posterior_array[i])\n",
    "        # find the postion of mean in x_points (closest position)\n",
    "        idx = np.argmin(np.abs(x_points - mean))\n",
    "        binary_decision_array[i, idx] = 1\n",
    "\n",
    "    return binary_decision_array\n",
    "\n",
    "binary_decision_array = calculate_binary_decision_array(x, posterior_array)\n",
    "with plt.xkcd():\n",
    "  plot_myarray(binary_decision_array,\n",
    "               'Chosen position $\\hat x$', 'Brain-encoded Stimulus $\\~ x$',\n",
    "               'Sample Binary Decision Array')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Section 5: Probabilities of encoded stimuli"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:57.614770Z",
     "iopub.status.busy": "2021-06-02T13:21:57.614295Z",
     "iopub.status.idle": "2021-06-02T13:21:57.705370Z",
     "shell.execute_reply": "2021-06-02T13:21:57.704892Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Video 5: Input array\n",
    "video = YouTubeVideo(id='C1d1n_Si83o', width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtube.com/watch?v=\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we as experimentalists can not know the encoding $\\tilde x$ of the stimulus $x$ that we do know, we had to compute the binary decision array for each possible encoding. \n",
    "\n",
    "First however, we need to calculate how likely each possible encoding is given the true stimulus. That is, we will now create a Gaussian centered around the true presented stimulus, with $\\sigma = 1$, and repeat that gaussian distribution across as a function of potentially encoded values $\\tilde x$. That is, we want to make a *column* gaussian centered around the true presented stimulus, and repeat this *column* Gaussian across all hypothetical stimulus values $x$.\n",
    "\n",
    "This, effectively encodes the distribution of the brain encoded stimulus (one single simulus, which we as experimentalists know) and enable us to link the true stimulus $x$, to potential encodings $\\tilde x$.\n",
    "\n",
    "**Suggestions**\n",
    "\n",
    "For this exercise, we will assume the true stimulus is presented at direction 2.5\n",
    "* Create a Gaussian likelihood with $\\mu = 2.5$ and $\\sigma = 1.0$\n",
    "* Make this the first column of your array and repeat that *column* to fill in the true presented stimulus input as a function of hypothetical stimulus locations.\n",
    "* Plot the array using the function `plot_myarray()` already pre-written and commented-out in your script"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Exercise 5: Generate an input as a function of hypothetical stimulus x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "code",
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:57.710828Z",
     "iopub.status.busy": "2021-06-02T13:21:57.709774Z",
     "iopub.status.idle": "2021-06-02T13:21:57.711395Z",
     "shell.execute_reply": "2021-06-02T13:21:57.711778Z"
    }
   },
   "outputs": [],
   "source": [
    "def generate_input_array(x_points, stim_array, posterior_array,\n",
    "                         mean=2.5, sigma=1.):\n",
    "\n",
    "    input_array = np.zeros_like(posterior_array)\n",
    "\n",
    "    ########################################################################\n",
    "    ## Insert your code here to:\n",
    "    ##      - Generate a gaussian centered on the true stimulus 2.5\n",
    "    ##        and sigma = 1. for each column\n",
    "    ## remove the raise below to test your function\n",
    "    raise NotImplementedError(\"You need to complete the function!\")\n",
    "    ########################################################################\n",
    "    for i in range(len(x_points)):\n",
    "        input_array[:, i] = ...\n",
    "\n",
    "    return input_array\n",
    "\n",
    "# Uncomment following lines, once the task is complete.\n",
    "# input_array = generate_input_array(x, hypothetical_stim, posterior_array)\n",
    "# plot_myarray(input_array,\n",
    "#              'Hypothetical Stimulus $x$', '$\\~x$',\n",
    "#              'Sample Distribution over Encodings:\\n $p(\\~x | x = 2.5)$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:57.716221Z",
     "iopub.status.busy": "2021-06-02T13:21:57.715751Z",
     "iopub.status.idle": "2021-06-02T13:21:58.332141Z",
     "shell.execute_reply": "2021-06-02T13:21:58.331706Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "def generate_input_array(x_points, stim_array, posterior_array,\n",
    "                         mean=2.5, sigma=1.):\n",
    "\n",
    "    input_array = np.zeros_like(posterior_array)\n",
    "\n",
    "    for i in range(len(x_points)):\n",
    "        input_array[:, i] = my_gaussian(stim_array, mean, sigma)\n",
    "\n",
    "    return input_array\n",
    "\n",
    "input_array = generate_input_array(x, hypothetical_stim, posterior_array)\n",
    "with plt.xkcd():\n",
    "  plot_myarray(input_array,\n",
    "               'Hypothetical Stimulus $x$', '$\\~x$',\n",
    "               'Sample Distribution over Encodings:\\n $p(\\~x | x = 2.5)$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Section 6: Normalization and expected estimate distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:58.336709Z",
     "iopub.status.busy": "2021-06-02T13:21:58.336229Z",
     "iopub.status.idle": "2021-06-02T13:21:58.485044Z",
     "shell.execute_reply": "2021-06-02T13:21:58.485449Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Video 6: Marginalization\n",
    "video = YouTubeVideo(id='5alwtNS4CGw', width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtube.com/watch?v=\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have a true stimulus $x$ and a way to link it to potential encodings, we will be able to calculate the distribution of encodings and ultimately estimates. To integrate over all possible hypothetical values of $\\tilde x$ we marginalize, that is, we first compute the dot-product from the true presented stimulus and our binary decision array and then sum over x. \n",
    "\n",
    "Mathematically, this means that we want to compute:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "    Marginalization Array = Input Array \\odot Binary Decision Array\n",
    "\\end{eqnarray}\n",
    "\n",
    "\\begin{eqnarray}\n",
    "    Marginal = \\int_{\\tilde x} Marginalization Array\n",
    "\\end{eqnarray}\n",
    "\n",
    "Since we are performing integration over discrete values using arrays for visualization purposes, the integration reduces to a simple sum over $\\tilde x$.\n",
    "\n",
    "**Suggestions**\n",
    "\n",
    "* For each row of the input and binary arrays, calculate product of the two and fill in the 2D marginal array.\n",
    "* Plot the result using the function `plot_myarray()` already pre-written and commented-out in your script\n",
    "* Calculate and plot the marginal over `x` using the code snippet commented out in your script\n",
    "   - Note how the limitations of numerical integration create artifacts on your marginal "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Exercise 6: Implement the marginalization matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "code",
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:58.489853Z",
     "iopub.status.busy": "2021-06-02T13:21:58.489369Z",
     "iopub.status.idle": "2021-06-02T13:21:58.492555Z",
     "shell.execute_reply": "2021-06-02T13:21:58.492110Z"
    }
   },
   "outputs": [],
   "source": [
    "def my_marginalization(input_array, binary_decision_array):\n",
    "\n",
    "    ############################################################################\n",
    "    ## Insert your code here to:\n",
    "    ##  - Compute 'marginalization_array' by multiplying pointwise the Binary\n",
    "    ##    decision array over hypothetical stimuli and the Input array\n",
    "    ##  - Compute 'marginal' from the 'marginalization_array' by summing over x\n",
    "    ##    (hint: use np.sum() and only marginalize along the columns)\n",
    "    ## remove the raise below to test your function\n",
    "    raise NotImplementedError(\"You need to complete the function!\")\n",
    "    ############################################################################\n",
    "\n",
    "    marginalization_array = ...\n",
    "    marginal = ... # note axis\n",
    "    marginal /= ... # normalize\n",
    "\n",
    "    return marginalization_array, marginal\n",
    "\n",
    "# Uncomment following lines, once the task is complete.\n",
    "# marginalization_array, marginal = my_marginalization(input_array, binary_decision_array)\n",
    "# plot_myarray(marginalization_array, 'estimated $\\hat x$', '$\\~x$', 'Marginalization array: $p(\\^x | \\~x)$')\n",
    "# plt.figure()\n",
    "# plt.plot(x, marginal)\n",
    "# plt.xlabel('$\\^x$')\n",
    "# plt.ylabel('probability')\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:58.497619Z",
     "iopub.status.busy": "2021-06-02T13:21:58.497144Z",
     "iopub.status.idle": "2021-06-02T13:21:59.195885Z",
     "shell.execute_reply": "2021-06-02T13:21:59.196292Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "def my_marginalization(input_array, binary_decision_array):\n",
    "\n",
    "    marginalization_array = input_array * binary_decision_array\n",
    "    marginal = np.sum(marginalization_array, axis=0)  # note axis\n",
    "    marginal /= marginal.sum()  # normalize\n",
    "\n",
    "    return marginalization_array, marginal\n",
    "\n",
    "marginalization_array, marginal = my_marginalization(input_array, binary_decision_array)\n",
    "with plt.xkcd():\n",
    "  plot_myarray(marginalization_array, 'estimated $\\hat x$', '$\\~x$', 'Marginalization array: $p(\\^x | \\~x)$')\n",
    "  plt.figure()\n",
    "  plt.plot(x, marginal)\n",
    "  plt.xlabel('$\\^x$')\n",
    "  plt.ylabel('probability')\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Generate some data\n",
    "\n",
    "We have seen how to calculate the posterior and marginalize to remove $\\tilde x$ and get $p(\\hat{x} \\mid x)$. Next, we will generate some artificial data for a single participant using the `generate_data()` function provided, and mixing parameter $p_{independent} = 0.1$. \n",
    "\n",
    "Our goal in the next exercise will be to recover that parameter. These parameter recovery experiments are a powerful method for planning and debugging Bayesian analyses--if you cannot recover the given parameters, something has gone wrong! Note that this value for $p_{independent}$ is not quite the same as our prior, which used $p_{independent} = 0.05.$ This lets us test out the complete model. \n",
    "\n",
    "Please run the code below to generate some synthetic data.  You do not need to edit anything, but check that the plot below matches what you would expect from the video. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:59.227732Z",
     "iopub.status.busy": "2021-06-02T13:21:59.225544Z",
     "iopub.status.idle": "2021-06-02T13:21:59.438751Z",
     "shell.execute_reply": "2021-06-02T13:21:59.438043Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title\n",
    "#@markdown #### Run the 'generate_data' function (this cell)\n",
    "def generate_data(x_stim, p_independent):\n",
    "  \"\"\"\n",
    "  DO NOT EDIT THIS FUNCTION !!!\n",
    "\n",
    "  Returns generated data using the mixture of Gaussian prior with mixture\n",
    "  parameter `p_independent`\n",
    "\n",
    "  Args :\n",
    "    x_stim (numpy array of floats) - x values at which stimuli are presented\n",
    "    p_independent (scalar) - mixture component for the Mixture of Gaussian prior\n",
    "\n",
    "  Returns:\n",
    "    (numpy array of floats): x_hat response of participant for each stimulus\n",
    "  \"\"\"\n",
    "  x = np.arange(-10,10,0.1)\n",
    "  x_hat = np.zeros_like(x_stim)\n",
    "\n",
    "  prior_mean = 0\n",
    "  prior_sigma1 = .5\n",
    "  prior_sigma2 = 3\n",
    "  prior1 = my_gaussian(x, prior_mean, prior_sigma1)\n",
    "  prior2 = my_gaussian(x, prior_mean, prior_sigma2)\n",
    "\n",
    "  prior_combined = (1-p_independent) * prior1 + (p_independent * prior2)\n",
    "  prior_combined = prior_combined / np.sum(prior_combined)\n",
    "\n",
    "  for i_stim in np.arange(x_stim.shape[0]):\n",
    "    likelihood_mean = x_stim[i_stim]\n",
    "    likelihood_sigma  = 1\n",
    "    likelihood = my_gaussian(x, likelihood_mean, likelihood_sigma)\n",
    "    likelihood = likelihood / np.sum(likelihood)\n",
    "\n",
    "    posterior = np.multiply(prior_combined, likelihood)\n",
    "    posterior = posterior / np.sum(posterior)\n",
    "\n",
    "    # Assumes participant takes posterior mean as 'action'\n",
    "    x_hat[i_stim] = np.sum(x * posterior)\n",
    "  return x_hat\n",
    "\n",
    "# Generate data for a single participant\n",
    "true_stim = np.array([-8, -4, -3, -2.5, -2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2,\n",
    "                      2.5, 3, 4, 8])\n",
    "behaviour = generate_data(true_stim, 0.10)\n",
    "\n",
    "plot_simulated_behavior(true_stim, behaviour)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "#Section 7: Model fitting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:59.443186Z",
     "iopub.status.busy": "2021-06-02T13:21:59.442668Z",
     "iopub.status.idle": "2021-06-02T13:21:59.596532Z",
     "shell.execute_reply": "2021-06-02T13:21:59.596900Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Video 7: Log likelihood\n",
    "video = YouTubeVideo(id='jbYauFpyZhs', width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtube.com/watch?v=\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have generated some data, we will attempt to recover the parameter $p_{independent}$ that was used to generate it.\n",
    "\n",
    "We have provided you with an incomplete function called `my_Bayes_model_mse()` that needs to be completed to perform the same computations you have performed in the previous exercises but over all the participant's trial, as opposed to a single trial.\n",
    "\n",
    "The likelihood has already been constructed; since it depends  only on the hypothetical stimuli, it will not change. However, we will have to implement the prior matrix, since it depends on $p_{independent}$. We will therefore have to recompute the posterior, input and the marginal in order to get $p(\\hat{x} \\mid x)$. \n",
    "\n",
    "Using $p(\\hat{x} \\mid x)$, we will then compute the negative log-likelihood for each trial and find the value of $p_{independent}$ that minimizes the negative log-likelihood (i.e. maximises the log-likelihood.  See the model fitting tutorial from W1D3 for a refresher).\n",
    "\n",
    "In this experiment, we assume that trials are independent from one another. This is a common assumption--and it's often even true! It allows us to define negative log-likelihood as:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "    -LL = - \\sum_i \\log p(\\hat{x}_i \\mid x_i)\n",
    "\\end{eqnarray}\n",
    "\n",
    "where $\\hat{x}_i$ is the participant's response for trial $i$, with presented stimulus $x_i$ \n",
    "\n",
    "* Complete the function `my_Bayes_model_mse`, we've already pre-completed the function to give you the prior, posterior, and input arrays on each trial\n",
    "* Compute the marginalization array as well as the marginal on each trial\n",
    "* Compute the negative log likelihood using the marginal and the participant's response\n",
    "* Using the code snippet commented out in your script to loop over possible values of $p_{independent}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Exercise 7: Fitting a model to generated data\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:59.604185Z",
     "iopub.status.busy": "2021-06-02T13:21:59.599207Z",
     "iopub.status.idle": "2021-06-02T13:21:59.606043Z",
     "shell.execute_reply": "2021-06-02T13:21:59.605643Z"
    }
   },
   "outputs": [],
   "source": [
    "def my_Bayes_model_mse(params):\n",
    "    \"\"\"\n",
    "    Function fits the Bayesian model from Tutorial 4\n",
    "\n",
    "    Args :\n",
    "        params (list of positive floats):  parameters used by the model\n",
    "        (params[0]  = posterior scaling)\n",
    "\n",
    "    Returns :\n",
    "        (scalar) negative log-likelihood :sum of log probabilities\n",
    "    \"\"\"\n",
    "\n",
    "    # Create the prior array\n",
    "    p_independent=params[0]\n",
    "    prior_array = calculate_prior_array(x,\n",
    "                                        hypothetical_stim,\n",
    "                                        p_independent,\n",
    "                                        prior_sigma_indep= 3.)\n",
    "\n",
    "    # Create posterior array\n",
    "    posterior_array = calculate_posterior_array(prior_array, likelihood_array)\n",
    "\n",
    "    # Create Binary decision array\n",
    "    binary_decision_array = calculate_binary_decision_array(x, posterior_array)\n",
    "\n",
    "    # we will use trial_ll (trial log likelihood) to register each trial\n",
    "    trial_ll = np.zeros_like(true_stim)\n",
    "\n",
    "    # Loop over stimuli\n",
    "    for i_stim in range(len(true_stim)):\n",
    "\n",
    "        # create the input array with true_stim as mean\n",
    "        input_array = np.zeros_like(posterior_array)\n",
    "        for i in range(len(x)):\n",
    "            input_array[:, i] = my_gaussian(hypothetical_stim, true_stim[i_stim], 1)\n",
    "            input_array[:, i] = input_array[:, i] / np.sum(input_array[:, i])\n",
    "\n",
    "        # calculate the marginalizations\n",
    "        marginalization_array, marginal = my_marginalization(input_array,\n",
    "                                                    binary_decision_array)\n",
    "\n",
    "        action = behaviour[i_stim]\n",
    "        idx = np.argmin(np.abs(x - action))\n",
    "\n",
    "        ########################################################################\n",
    "        ## Insert your code here to:\n",
    "        ##      - Compute the log likelihood of the participant\n",
    "        ## remove the raise below to test your function\n",
    "        raise NotImplementedError(\"You need to complete the function!\")\n",
    "        ########################################################################\n",
    "\n",
    "        # Get the marginal likelihood corresponding to the action\n",
    "        marginal_nonzero = ... + np.finfo(float).eps  # avoid log(0)\n",
    "        trial_ll[i_stim] = np.log(marginal_nonzero)\n",
    "\n",
    "    neg_ll = - trial_ll.sum()\n",
    "\n",
    "    return neg_ll\n",
    "\n",
    "# Uncomment following lines, once the task is complete.\n",
    "# plot_my_bayes_model(my_Bayes_model_mse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T13:21:59.613324Z",
     "iopub.status.busy": "2021-06-02T13:21:59.612795Z",
     "iopub.status.idle": "2021-06-02T13:22:13.129015Z",
     "shell.execute_reply": "2021-06-02T13:22:13.128580Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "def my_Bayes_model_mse(params):\n",
    "    \"\"\"\n",
    "    Function fits the Bayesian model from Tutorial 4\n",
    "\n",
    "    Args :\n",
    "        params (list of positive floats):  parameters used by the model\n",
    "        (params[0]  = posterior scaling)\n",
    "\n",
    "    Returns :\n",
    "        (scalar) negative log-likelihood :sum of log probabilities\n",
    "    \"\"\"\n",
    "\n",
    "    # Create the prior array\n",
    "    p_independent=params[0]\n",
    "    prior_array = calculate_prior_array(x,\n",
    "                                        hypothetical_stim,\n",
    "                                        p_independent,\n",
    "                                        prior_sigma_indep= 3.)\n",
    "\n",
    "    # Create posterior array\n",
    "    posterior_array = calculate_posterior_array(prior_array, likelihood_array)\n",
    "\n",
    "    # Create Binary decision array\n",
    "    binary_decision_array = calculate_binary_decision_array(x, posterior_array)\n",
    "\n",
    "    # we will use trial_ll (trial log likelihood) to register each trial\n",
    "    trial_ll = np.zeros_like(true_stim)\n",
    "\n",
    "    # Loop over stimuli\n",
    "    for i_stim in range(len(true_stim)):\n",
    "\n",
    "        # create the input array with true_stim as mean\n",
    "        input_array = np.zeros_like(posterior_array)\n",
    "        for i in range(len(x)):\n",
    "            input_array[:, i] = my_gaussian(hypothetical_stim, true_stim[i_stim], 1)\n",
    "            input_array[:, i] = input_array[:, i] / np.sum(input_array[:, i])\n",
    "\n",
    "        # calculate the marginalizations\n",
    "        marginalization_array, marginal = my_marginalization(input_array,\n",
    "                                                    binary_decision_array)\n",
    "\n",
    "        action = behaviour[i_stim]\n",
    "        idx = np.argmin(np.abs(x - action))\n",
    "\n",
    "        # Get the marginal likelihood corresponding to the action\n",
    "        marginal_nonzero = marginal[idx] + np.finfo(float).eps  # avoid log(0)\n",
    "        trial_ll[i_stim] = np.log(marginal_nonzero)\n",
    "\n",
    "    neg_ll = - trial_ll.sum()\n",
    "\n",
    "    return neg_ll\n",
    "\n",
    "with plt.xkcd():\n",
    "  plot_my_bayes_model(my_Bayes_model_mse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Section 8: Summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T13:22:13.133803Z",
     "iopub.status.busy": "2021-06-02T13:22:13.133273Z",
     "iopub.status.idle": "2021-06-02T13:22:13.263917Z",
     "shell.execute_reply": "2021-06-02T13:22:13.263508Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Video 8: Outro\n",
    "video = YouTubeVideo(id='F5JfqJonz20', width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtube.com/watch?v=\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Congratuations! You found $p_{independent}$, the parameter that describes how much weight subjects assign to the same-cause vs. independent-cause origins of a sound. In the preceeding notebooks, we went through the entire Bayesian analysis pipeline:\n",
    "\n",
    "*   developing a model\n",
    "*   simulating data, and\n",
    "*   using Bayes' Rule and marginalization to recover a hidden parameter from the data\n",
    "\n",
    "This example was simple, but the same princples can be used to analyze datasets with many hidden variables and complex priors and likelihoods. Bayes' Rule will also play a cruical role in many of the other techniques you will see later this week. \n",
    "\n",
    "---\n",
    "\n",
    "If you're still intrigued as to why we decided to use the mean of the posterior as a decision rule for a response $\\hat{x}$, we have an extra Bonus Tutorial 4 which goes through the most common decision rules and how these rules correspond to minimizing different cost functions."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "W3D1_Tutorial3",
   "provenance": [],
   "toc_visible": true
  },
  "kernel": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
